A WENO Algorithm for the Radiative 
Transfer and Ionized Sphere at Reionization 



Jing-Mei Qiu^, Chi- Wang Shu^, Long-Long Feng^''^, 

Li-ZhiFang^ 

\o ■ 

f^ ' ^Division of Applied Mathematics, Brown University, Providence, RI 02912, USA 

^^ ' Purple Mountain Observatory, Nanjing, 210008, P.R. China 

^ ; 

CIh, '^National Astronomical Observatories, Chinese Academy of Science, Chao-Yang 

<^ ; District, Beijing 100012, P.R. China 

CN \ Department of Physics, University of Arizona, Tucson, AZ 85721, USA 



> 

OO ; Abstract 

(N 

■^ ' We show that the algorithm based on the weighted essentially nonoscillatory 

^^ ■ (WENO) scheme with anti-diffusive flux corrections can be used as a solver of the 

^) \ radiative transfer equations. This algorithm is highly stable and robust for solving 

"^ I problems with both discontinuities and smooth solution structures. We test this code 

Q^' with the ionized sphere around point sources. It shows that the WENO scheme can 



o 



reveal the discontinuity of the radiative or ionizing fronts as well as the evolution of 



H \ photon frequency spectrum with high accuracy on coarse meshes and for a very wide 

^ ■ parameter space. This method would be useful to study the details of the ionized 

patch given by individual source in the epoch of reionization. We demonstrate this 

method by calculating the evolution of the ionized sphere around point sources in 

/\ ' physical and frequency spaces. It shows that the profile of the fraction of neutral 

j^ ■ hydrogen and the ionized radius are sensitively dependent on the intensity of the 

source. 
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1 Introduction 



In the study of cosmic large scale structures formation, one generally focuses 
on the inhomogeneity of the density and velocity fields of dark matter and 



Preprint submitted to Elsevier Science 5 February 2008 



baryon gas, but assumes that the distribution of background radiation is uni- 
form. This would be a reasonable approximation when the universe is trans- 
parent. Therefore, the effect of radiative transfer is generally ignored in the 
study of large scale structures on low redshift. However, the fluctuations in 
the radiation field is essential in understanding the early universe, especially 
in the epoches around reionization. Radiative transfer problems are no longer 
to be overlooked. Many numerical solvers for the radiative transfer equation 
have been proposed (Abel et al. 1999, Ciardi et al. 2001, Gnedin & Abel 2001, 
Sokasian et al. 2001, Razoumov et al. 2002, Gen 2002, Maselh et al. 2003). 

In this paper, we will introduce the algorithm of the radiative transfer equa- 
tions based on the finite difference weighted essentially non-oscillatory (WENO) 
scheme (Jiang & Shu 1996) modified with anti-diffusive flux corrections (Xu 
& Shu 2005). It has been shown that the Boltzmann equations can be solved 
by this WENO method with reasonable computational speed at least for one 
or two spatial dimensions and two or three phase space dimensions plus time 
(Garrillo et al. 2003, 2006). The WENO code is effective to capture disconti- 
nuities as well as to resolve complicated flow features. Moreover, it is also sig- 
nificantly superior over piecewise smooth solutions containing discontinuities. 
It has advantages over the Monte-Garlo codes in eliminating fluctuations, in 
obtaining accurate probability density functions, and in accurately resolving 
time transients of the system. Gonsidering the equation of radiative trans- 
fer basically is the same as the Boltzmann equation, the WENO algorithm 
for the radiative transfer is straightforward. Moreover, a hybrid algorithm of 
hydrodynamic/N-body simulation based on WENO scheme has been success- 
fully developed (Feng et al. 2004). This result also motivated us to extend 
the WENO scheme to radiative transfer problems. It would be a necessary 
preparation to develop solver of hydro dynamic /radiative transfer problems. 

We will demonstrate the WENO algorithm by applying the WENO code to 
the Stromgren sphere. The structure and evolution of an individual Stromgren 
sphere are important for the understanding of the reionization process. The 
ionization of baryon gas surrounding high redshift objects, such as quasars, 
stars of first generation, generally is modelled by the Stromgren sphere (Madau 
& Rees, 2000; Gen & Haiman 2000; White et al. 2003; Stuart et al. 2004; Yu & 
Lu 2005). The Stromgren sphere usually is described as a fully ionized region, 
sharply divided from the neutral hydrogen region on the outside of the sphere. 
Therefore, we should require the algorithm to be able to 1.) capture the sharp 
boundary between the ionized and neutral regions, and their time-dependence; 
2.) accurately calculate the small amount of neutral hydrogen in the sphere 
and ionized hydrogen outside the sphere. These two points can be realized by 
the the WENO numerical scheme, as it is effective for problems containing 
both strong discontinuities and smooth flow structures. 

The paper is organized as follows. Section 2 describes the WENO numerical 



scheme. Section 3 gives the numerical solutions of the Stromgren sphere. The 
discussion and conclusion are presented in §4. The derivation of the equation 
of the radiative transfer is given in the Appendix. 



2 Numerical solver of the WENO scheme 



The radiative transfer equation in an expanding universe is (see the Appendix) 



^^ ^ ^''^-j]+^iHJ) = -iK + 3H)J + S, 



d (ct) dx^ \ '^ / '^^ 

where J{t,:si,i',ni) is the specific intensity, a the cosmic factor, H = a/a, v 
the frequency of photon, lo = Inl/z/, and n, a unit vector in the direction of 
photon propagation, /c^ and S are, respectively, the absorption and sources of 
photons. 

To demonstrate the WENO method, we consider the simplest case. We ignore 
the expansion of the universe. This approximation is correct if the time scale 
of evolution of photon distribution is much less than the time scale of cosmic 
expansion. Moreover, the physical and phase spaces are assumed to be 1- 
dimensional with coordinate r. In this case, the radiative transfer equation 
is 

The computational domain is discretized into a tensor product mesh. The 
mesh is taken to be uniform in the r-direction, and to be smooth non-uniform 
in the zz-direction, i.e. 



iAr; i = 0,1,2, ..., Nr, 



V 



3 



:2«- J =0,1, 2,..., AT,, (3) 



'max 



where Ar = r^s.x/Nr is the mesh size in the r-direction, and ^j = jA^, A^ 
loggZ^max/A^^ is the transformed mesh size in the z/ direction, r^ax and u^ 
are, respectively, the sizes of the numerical domain, which are adjusted in the 
numerical experiment such that J{t,r,v) ~ 0, for r > rmax and all t and v, 
and J{t, r, u) ~ 0, for z/ > z/max and all t and r. 

The approximations to the point values of the solution J(t"',ri,i'j), denoted 
by Jf-, are obtained with an approximation to the spacial derivatives using 
the 5th order finite difference WENO method (Jiang & Shu 1996) with anti- 
diffusive flux corrections (Xu & Shu 2005). 



Approximation to the derivatives: 

To calculate dJ /dr, the variable v is fixed and the approximation is per- 
formed along the r-line 

(9 1 ^ ^ 

where the numerical fiux h1_^_^,2 is obtained with the procedure given be- 
low. We can use the upwind fiuxes without fiux splitting in the fifth order 
WENO approximation because the wind direction is fixed (positive). To 
obtain the sharp resolution of the contact discontinuities, the anti-diffusive 
fiux corrections are used in our code. 
First, we denote 

h, = J(t", n, Uj),i = -2, -1, ..., Nr + 2 

where n and j are fixed. The numerical fiux from the regular WENO pro- 
cedure is obtained by 

where hj^l/2 are the three third order fiuxes on three different stencils given 

by 

'^i+l/2 — o '^J-2 — « '^*-l + ~fi~ ' 



,(2) _ 1. 5. 1 

/2- 

1, 5, 1 



^i+l/2 ~ ~R^*~1 ~^ R^* ~^ o^J+1' 



and the nonlinear weights uj^ are given by 



UJm — TvH ^? ^l 



with the linear weights 7^ given by 



71 ^Q, 72 g, 73 ^Q, 

and the smoothness indicators /3/ given by 

A = — (/i^-2 - 2/ii_i + /i,)' + - (/i,_2 - 4/ii_i + 3hif 
13 2 1 2 

/52 = — (/ij-l - 2/lj -h /li+l) + - (/li_l - /li+l) 

13 2 1 2 



e is a parameter to avoid the denominator to become and is taken as 
e = 10"^ times the maximum magnitude of the initial condition J in the 
computation of this paper. The reconstruction of the finite difference WENO 
flux on the downwind side hf^-^,^ is obtained in a mirror symmetric fashion 

with respect to Xi+1/2 as that for /ij^i^- 

The anti-diffusive flux corrections are based on the fluxes obtained from 
the regular WENO procedure. It is given by 

ht+i/2 = Ki/2+ (4) 



hi — /ij_l , , „ , _ ; ^ 



(/.^minmod I + h^^^^^ " ^*+i/2> K+1/2 " ^^+1/2 I > 

where 77 = At/Ar is the cfl number and the minmod function is defined as 



minmod(a, b) 



0, if ab<0 

a, if ah > 0,\ a \<\ b \ (5) 

b, if a6 > 0, 1 6 |<| a I . 



(pi of eq.(4) is the discontinuity indicator between and 1, defined as 

Pi + li 
where 

I |2 

Pi = { \ ), li = , ai = {\ hi^i - hi \ +0 , 

(Xi-\ ai+2 Oii 

with (^ being a small positive number taken as 10^^ in our computation. 
■Umax and -Umin a^G the maximum and minimum values of hi for all grid 
points. With the definition above, we will have < 0j < 1. 0j = O(Ar^) 
in the smooth regions and 0j is close to 1 near strong discontinuities. The 
purpose of the ant i- diffusive flux corrections is to improve the resolution 
of contact discontinuities without sacrificing accuracy and stability of the 
original WENO scheme. 

Approximation to the evolution with time 

The time-derivative d/dt in each step At is calculated by the third order 
TVD Runge-Kutta method as 

J« = J" + AtL(J",r) 

j(2) _ jn ^ lAtL(J") + iAtL(/i)) 

jn+i ^jn^ -^tLi.r) + -AtL( J(i)) + -AtLU^^^) 
6 6 3 



where L is the approximation to the spatial derivatives and the source terms: 



L{J)^-^J-k,J + S 
or 

The Runge-Kutta method needs to be modified considering the modifi- 
cation on the anti-diffusive flux /" by 

J« = J" + AtL(J",r) 

j(2) _ jn ^ -^tL'ir) + -AtL( J(^)) 
4 4 

jn+l ^ jn ^ \^tL"^jn) + -ML{J^^'^) + -AtL{J^^^) 

6 6 3 

where the operator L is deflned by the anti-diffusive flux /i" given by eq.(4), 
and the operator L' is deflned by the modifled anti-diffusive flux h as 



K 



^j+i/2 + niinmod 






J+l/2 



+"'j-l/2 "-4+1/2' "'j+1/2 "-4+1/2^' 

?/ 6c> 0, 1 6 |<| c I, 
otherwise 



'^i+l/2; 



(6) 



and L" is deflned by the modifled anti-diffusive flux h^, 

'6{hi - hi^i) 



h. 



^i+i/2 + minmod 



T] 



i+l/2 



"'"'^j-1/2 '''i+l/2-' '''i+l/2 '''i+l/2 ) 1 

if 6c> 0, 1 6 |<| c 
otherwise 



"'i+l/2^ 



(7) 



Here b = {hi - hi.i)/r] + h._y^ - ^^^^^ c = hj^^y^ - \_^y^. 

We have given the details of the WENO algorithm with anti-diffusive flux 
corrections only for the one dimensional case to save space. The flnite dif- 
ference WENO scheme is ideally suited for multi-dimensional calculations, 
as derivatives in each direction can be approximated in an one dimensional 
setting by flxing the other variables, while still maintaining high order accu- 
racy and stability. We refer to (Carrillo et al. 2003, 2006) for more details of 
WENO approximations to Boltzmann equations. Our radiative transfer equa- 
tions [Eq.(l)] are of the same form as the Boltzmann equations in (Carrillo et 
al. 2003, 2006) and hence the algorithms developed there can be applied here 
without difficulty. 



3 Ionized sphere 

3.1 Ionized source in a uniform medium 



As a test, we consider a point plioton source located at tlie center x = of 
a uniformly distributed medium. Assuming the source is monochromatic with 
frequency vq, and all photons are emitted along the radial direction, we have 



S{t, X, z/, n) = /(t, x)5(z/ - Vo)5{n - e^). 



and 



, ElV, at the center x = 
/(t,x) = <{ (9) 

0, otherwise. 



where E is the total energy of photons emitted from the sources per unit time. 
When l^ -^ 0, /(t,x) -^ E6{x.). 

In this case, Eq.(l) could be written in the spherical coordinates as 

dJ .19/2, 



+ ^7- ir'n'J) = -kJ, r ^ (10) 



d (ct) r^ dr 



where the absorption coefficient k is given by A; = na, n being the number 
density of particle of the medium, and a the absorption cross-section. In this 
problem, both n and a are assumed to be constants. 

We define r'^J = J'{t,r)5{u — Uo)5{yi — ey.). Then eq.(lO) becomes 

dJ' d 

+ ^f = -kJ', r^O. (11) 



d (ct) dr 



For simplicity, we drop the prime, and use J(t, r) for J'(t, r) below. It will not 
cause confusion between the J in eqs.(ll) and (10), as the former is a function 
of t, r, while the later is a function of t, r, u and rij. The source term gives a 
boundary condition at r = as 

\im4:nJ{t,r) =E (12) 
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Fig. 1. The curves of J{t,r) given by the numerical solution at Nr = 200 (square) 
and the exact solution (solid line). J(t=0, r)=^. The time is taken to be t = 1.0 
(top left); t = 2.0 (top right); t = 3.0 (bottom left) and t = 4.0 (bottom right). The 
small windows in the bottom figures are the zoom-in version of the corresponding 
regions in the circles. 

Assuming the source starts to emit photon at t = 0, the initial condition is 
then 



J(t = 0,r) 




(13) 



Subject to the conditions (12) and (13), the exact solution of eq.(ll) is 

E 



J(t,r) = 9{t-r)—e 



—kr 



(14) 



where d{x) is a step function: 9{x) = 1, if a; > 0; and 9{x) = 0, if a: < 0. The 
step function is from the radiative front, which propagates with the speed of 
light r = t. 

Eq.(ll) has the same shape as Eq.(2). Similarly, we can apply the WENO 
scheme to solve Eq.(ll). The boundary condition eq.(12) makes us to take the 



inflow condition at r = as 

J{f,n) = E/ATi, for- ^ = 0,-1, -2. 
The boundary condition at r = rmax is: 

JNr+i,j = JNr-i,j, for i = 0, 1, 2 

Figure 1 plots both the numerical solution and exact solution at time t = 1.0, 
2.0, 3.0, 4.0. Clearly, the numerical result displays an excellent agreement with 
the exact solution even when the number Nr is only 200. The discontinuity at 
the radiative front r = t is also well reproduced in the WENO scheme. 

3.2 The profile of an ionized sphere 



We now consider a realistic problem: the ionization of a uniformly distributed 
hydrogen gas with number density n by a point UV photon source located at 
r = 0. This problem is the well known Stromgren sphere if the proflle of the 
ionized sphere is approximated by two sharply divided regions as 

, if r < i?, 

/Hi(r) ^ (15) 

1 if r > i?. 



where /hi = rim/n, nHi(t,x) is the number density of neutral hydrogen, HI. 
Eq.(15) means that within the sphere of radius Rg around the source, hydrogen 
is fully ionized, while outside the sphere hydrogen atoms remain neutral. Rs is 
called the ionized radius of the Stromgren sphere, which can be determined by 
the balance between the rate of recombinations and the emission N of ionizing 
photons (Stromgren 1939, Spitzer 1978, Osterbrook 1989) 



1/3 

Rs = ( ._r' „. ) , (16) 



/ 3N 
I 47rQ;Hii^^ 



where anii is the hydrogen recombination coefficient. 

The sharp proflle eq.(15), and then eq.(16), are reasonable if the mean free 
path of photon A ~ l/o'o^ is much less than Rg. For strong sources, like 
quasars A^ ~ 10^^ s^^, Rg is much larger than A. However, for weak sources, 
like population HI stars A^ < 10^°, A is comparable, or even larger than Rg. 
In this case, the proflle eq.(15) is no longer valid. Even for strong sources, 
we need to study the small deviation of the proflle of ionized sphere from 



eq.(15). The profile should be found by solving the radiative transfer equation 
with proper boundary and initial conditions. Moreover, the profile eq.(15) 
describes the final state of the ionized sphere around a strong source. To 
study the reionization history of the universe, we may need the information of 
the formation of the ionized sphere i.e. the time-dependence of fuiiit, r). This 
also requires to solve the radiative transfer equation. 

To calculate the profile of the ionized sphere, we can still use eq.(ll), boundary 
condition eq.(12), and initial condition eq.(13). The absorption coefficient now 
is given by 

/f^i^o = ^('^o)"'Hi(i,x) (17) 



where o"(z/o) is the absorption cross section at the ionization frequency uq. We 
have o-Q = o-(z/o) = 6.3 x 10"^^ cm^. 

The number density of neutral hydrogen nHi(t, x) is determined by the ion- 
ization equilibrium equation 



dfm _ 

—— — CtHII^e/HII — -L 7HIJHI — -L cHI^e /HI- 



where fui(t,r) = riui/n, f}iu(t,r) = nuii/n. We will assume electron density 
He = hhii- This means that the electrons from ionized helium are ignored. The 
photoionization is given by 

1 J(t,r) , , 



The parameters anii and FeHi are the recombination coefficient and collision 
ionization, respectively, which can be found in e.g. Theuns et al (1998). 

Let us rescale the variables by t' = co'{i'o)nt, r' = cr{i'o)nr. It means that t' and 
r' are, respectively, the time and distance in the units of mean free flight time 
and mean free path, l/(To(t'o)^) of photon huQ in neutral hydrogen gas with 
density n. In the ACDM model, n = 1.88 x 10~^(1 + 2;^)^ cm~^, where Zr is the 
redshift of reionization, the unit oft' is 0.89 x 10^(1 + 2;^)"^ years, and the unit 
of r' is 0.27(1 + 2;^)"^ Mpc. We also rescale the intensity by J" = J{a'^n/chi'o). 
Eqs.(ll) can then be rewritten as 

(9 T" d T" 

10 



The boundary condition (12) and initial condition (13) are, respectively 

j"(t,r = 0) = J^ (21) 



J"(t = 0,r) 



0, 
7" 



if r > 
if r = 



(22) 



With this rescaling, the photoionization term of eq.(19) is r^Hi/"^ = ( J"/T'^)(cno'o) 
The intensity of source is iV = 47rJ(0)//iz/ = 5.05 x 10^Vo'(l + Zr)"^ s^^ 



3.2.1 Strong sources 

For strong sources, Rg is much larger than the mean free path. We can replace 
the rate equation (18) by ionization equilibrium equation dfm/dt = 0, or 
ttHii/m ~ (X'y,ui/n + 2aHii)/H/ + ctHii = 0, where the small term Fem is also 
dropped. A typical numerical result of the fuiit, r) profile is shown in Figure 2, 
in which Jq = 4.2, or N = 2.16 x 10^^(1 + Zr)~^. The calculation is performed 
with Nr = 1000. 
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Fig. 2. /Hi(t, r) vs. r at time t = 30.0 (top left), 60.0 (top right), 90.0 (bottom left), 
140.0 (bottom right) for a strong source Jq = 4.2, or iV = 2.16 x 10^^(1 + Zr)^"^ 
■sT^ . The parameter N^r of the numerical calculation is taken to be 1000. The ionized 
radius i?s=122.04. 

Figure 2 shows that the code can well reveal the jump at radiative front r = t, 
which propagates with the speed of light. The evolution of the ionized range 



11 



can approximately be described as 



fm{t,r) 



9{r-t), t < ts 
9{r-ts),t>ts 



(23) 



When t is small, the increase of ionized sphere is following the radiative front 
r = t. The increase will be halted, when the ionization equilibrium is totally 
established, and the ionized sphere becomes time-independent. For the case 
of J"(0) = 4.2, the final state arrives at tg ^ 1.22 x 10^. Since the bottom 
right panel of Figure 2 is for t = 140 > t^, it should be the final state of the 
ionized sphere. The profile of the final state can be approximately described 
by the Stromgren sphere profile eq.(15), and ionized radius Rg is about 122 
mean free path. 
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Fig. 3. /Hi(t,r) vs. r at time t = 400.0 (top left), 800.0 (top right), 1000.0 
(bottom left), 1100.0 (bottom right) for a strong source Jq" = 4.2 x 10^ or 
iV = 2.6 X 10^^(1 + ^^r)^^ s^-*^. The parameter A^^ of the numerical calculation 
is taken to be 200. The ionized radius i?s=1080.0. 

If we define an effective ionized radius rHii(t) by 



r=t 



-r 



nil 



{t)= J fmi{r,tydr, 



(24) 



we have rHii(t) ^ t when t < tg- The physical meaning of rnn is clear. It gives 
an equivalent sphere, within which hydrogen is fully ionized. The ionized radius 
is essentially the same as the radiative front. Therefore, for strong sources, the 
radiative front and ionized radius are the same in all time t < tg- 



12 



From Figure 2 we can also see that the neutral hydrogen fraction within the 
ionized sphere r < tg actually is not exactly equal to zero. Although the profile 
(15) may still be reasonable, /hi can be as large as 10~^ near the ionized radius. 
We also calculated the profile with very strong source like Jq = 4.2 x 10^, which 
is shown in Figure 3. In this case the source intensity Jq is larger than that of 
Figure 2 by a factor of 10^, and therefore, the radius of Rg of Figure 3 is larger 
than that of 2 by a factor of 10. There is also neutral hydrogen left behind 
the ionization front. These small amounts of HI are stable, i.e. independent 
of the parameter A^^- The neutral hydrogen remained in the ionized sphere is 
not always negligible. 



3.2.2 Weak Sources 

If the Stromgren sphere radius Rs of a source is comparable or less than the 
mean free path, it is a weak sources. In this case, we should use eq.(18) to 
describe the ionization evolution. In the range of r less than the mean free 
path, i.e. r' <^ 1, the photonionization Fhi is large. We can keep only the 
photoionization term Fhi/hi on the r.h.s. of eq.(18). 




Fig. 4. /Hii(t,r) vs. r at time t = 1.0 (top left), 2.0 (top right), 3.0 (bottom left), 
4.0 (bottom right) for a weak source J^' = 0.001 or iV = 5.05 x 10^*^(1 + Zr)~^. The 
parameter iV^ of the numerical calculation is taken to be 200. 

A solution numerical solution of the profile /Hii(i, r) = 1 — fui(t, r) at t=l, 2, 
3, and 4 are shown in Figure 4, in which the source intensity is Jg = 0.001 or 
N = 5.05 X 10'^^(1 + Zr)^^ s"^ The calculation is performed with Nr = 200. 

We can see from Fig. 4 that the profile is very different from strong source. 
The spatial range of fm(t, r) significantly less than 1 (or fmiit, r) significantly 
larger than zero) is much less than t, because the total number of photons 



13 



emitted from sources within time t is less than the number of atoms within 
radius r = t. /Hi(t, r) is gradually increasing with r. Although the non-zero 
range of the ionized fraction /hh = 1 — fmii-, f^) increases with time, the range 
/hi(^)^) > 10~^ is always small. It implies that the ionized sphere probably is 
not transparent to Lja photons. 

Figure 4 shows that the effective ionized radius rHii(t) [eq.(24)] for weak 
sources is rHii(t) < t, or rHii(t)/t < 1. The growth of the ionized range is 
much slower than the radiative front. With the increase of source intensity, 
^Hii(^)/^ will gradually approach to 1, and then, the solution will transfer to 
the case of strong sources. 



3.3 Evolution of the frequency spectrum 



Let us consider a point source emitting photons with a power law frequency 
spectrum. The source function S is then 

(E{iy)/V)6{n- Br), at center 
5(t, X, z/, n) = <( (25) 

0, otherwise. 



where E^u) = E[v / vq)~'^ , a is the index of the power law. Therefore, different 
from the problems of §3.1 and 3.2, we should consider the variable v of the 
frequency space, i.e. J is a function of t, r, and v. The equation of J is still 
the same as eq.(ll), but the absorption coefficient is //-dependent 

ky = o-(z/)nHi(t,x), 

where the cross section cr[i>) = 6.3x 10^^^(z/o/^)^ cm^. The boundary condition 
now should be 

limJ(t,r,z/) = Jo(z//z/o)"" (26) 

r— >0 



and the initial condition is 

J{t = 0, r, u) = Jo(z//z/o)""^(r = 0) (27) 

We still rescale the variables by t' = cna'{i'o)t, r' = no'(z/o)r, u' = v/vq and 
J" = Jna'^{i'o)/ch. Then eq. (19) is updated as follow, 

(9 7" dJ" / 1 \ ^ 

^ ^ fmJ" (28) 



dt' dr' \ V / 
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Fig. 5. A numerical result of the intensity J"(t, r, u) for Jq = 0.001 at time t = 1.0 
(top left), 2.0 (top right), 3.0 (bottom left), and 4.0 (bottom right). The mesh is 
taken to be iV^=100, iV,,=100, and fmax=10'^. 

If we ignore the ionizing heating, /hi or nHi(t, x) is still determined by eq.(18), 
but r^Hi(^,x) should be given by 



oo 

r^Hi(t,x) = — j dv 



J{t,r,u) 



hv 



-aiu . 



(29) 



1^0 



This integration could be worked out with a quadrature formula accurate of 
order 4 



^ oo 

f{x)dx = AxY,WjfU^x) + 0{Ax^ 



VO 



.?=.70 



(30) 



where z/q = joAx, and the weights Wj are given by 



w 



JO 



, 5 "^JO 



W,o + l 



7 23 

6' ^^«+^ = ^' 



and 

Wja+j = 1> for j > 2. 

By having a non-uniform mesh in the //-direction, one cannot use this quadra- 
ture formula directly. However, since a uniform mesh size is used on ^, with 
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^3 



2^\ we can perform the numerical integration with respect to C, using 



(29). 




Fig. 6. A numerical solution of J vs. log v with A^^ = 100, Ni, = 100 (square) and the 
reference solution (solid) obtained with A^^ = 1000, A^;/=1000, at time t = 4.0 and 
r=0.8 (top left), r=2.0 (top right), r=2.8(bottom left) and r=3.6 (bottom right). 



^ and a = 1 is shown in Figure 
,=100 and t'max=10®. It displays 



A numerical result of J"{t,r,iy) for Jq = 10" 

5, which is obtained by the mesh N^=100, Ni 

the evolution of the frequency-dependence of J". We notice that for the high 

frequency u > 50, the intensity J" is almost independent of t and r, while the 

frequency spectrum at z/ < 10 is significantly dependent on both t and r. 
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Fig. 7. a vs. log u of the numerical solution with A^^ = 100, N^ = 100 at time t=4.0 
and r=0.8 (dash dot dot line), r=2.0 (dash line), r=2.8 (dash dot line) and r =3.6 
(solid line). 

The evolution of the frequency spectrum can be seen more clearly in Figure 
6, which shows J vs. u at time t = 4.0 and position r = 0.8, 2.0, 2.8 and 
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Fig. 8. a vs. r with samples of iV^ = 100, N^ = 100 at time t = 1.0 (left) and 4.0 
(right) for wave bands u = 1.32 (solid line), i^ = 1.74 (dash line), u = 2.29 (dash 
dot line), v = 3.02 (dot line), v = 3.98 (long dash line) and v = 63.10 (dash dot dot 
line). 

3.6. The results of Figure 6 given by A^^. = 100 and A^j, = 100 are the same as 
that obtained with A^^ = 1000 and N,^ = 1000. Therefore, the algorithm is also 
stable and highly accurate in the frequency space. From Figure 6, we can see a 
significant r-dependence of the spectrum. At small r, the spectrum essentially 
is of power law, while at r = 3.6 it is similar to a spectrum of self-absorption, 
i.e. it is low at z/ ~ 1, and shows a peak at z/ ~ 2. 

Generally, self- absorption leads to the lack of photons with z/ ~ 1, and there- 
fore, to the hardening of the frequency spectrum of photons. One can measure 
the hardening by the index of power law defined by 



a 



din J 
91nz/ 



(31) 



Figure 7 plots a vs. u at different positions r and shows that a becomes smaller 
at z/ < 10. Figure 8 gives a vs. r for several z/. It shows that a becomes smaller 
at larger r. 



4 Concluding remarks 



We described a numerical solver for radiative transfer problems based on the 
WENO scheme modified with ant i- diffusive flux corrections. It has high order 
of accuracy and good convergence, and is also highly stable and robust for 
solving problems with both discontinuities and smooth solution structures. 
Using the Stromgren sphere or ionized sphere as numerical tests, we showed 
that this code is able to resolve the sharp ionized front under a wide parameter 
range, such as intensity of the source A^ varying from 10^^ — 10^^(1 + Zr)~^ 
s~^, which covers various sources responsible for the early reionization in the 
universe. 
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Since the WENO scheme needs more floating point operations per cell than 
those of the PPM and TVD schemes, it leads to twice or more loss of com- 
putational speed. Nevertheless, it has been already successfully applied to 
kinetic equations of distribution function in phase space with one or two spa- 
tial dimensions and two or three phase space dimensions. It can provide the 
solution of evolution in both the physical space as well as the frequency space. 
Therefore, the WENO scheme would be useful to study the early stage of 
reionization, in which ionized region is generally an isolated patch given by 
individual source with various configuration. 



We demonstrate this algorithm with the problems of spherical ionized region 
exposed to a point source. It is shown that the radial profile of the fraction 
of ionized hydrogen is sensitively dependent on the intensity of the sources. 
Although for strong sources the Stromgren sphere provides a good approxi- 
mation to the ionized sphere, the profile of the sphere actually can not always 
be approximated as fully ionized spheres sharply separated with neutral hy- 
drogen region. Even for strong sources, the neutral hydrogen remained in the 
Stromgren ionized sphere may still be significant and important. There are 1-D 
radiative transfer codes, such as CLOUDY94 (e.g. Maselli et al. 2003), which 
are capable of capturing the sharp boundary of the Stromgren sphere. How- 
ever, they yield large errors of the neutral hydrogen remained in the sphere. 
Our result shows that the WENO codes can effectively handle the uncertainty 
around discontinuities. 



We studied the time-dependence of the ionized radius rHii(t). We found that 
for strong sources, rHii(t) — t when rnn < Rg- This is similar to the result 
of Madau & Rees 2000. For weak sources, we have rHii(t) < t. However, 
the evolution of rHii(t) can not be fitted by rHii(t) oc t^^^, which will yields 
unphysical result dri{ii(t)/dt > 1 at small t (e.g. Cen & Haiman 2000; Yu & 
Lu 2005). Moreover, the WENO algorithm provides also the solution of the 
evolution of the photon's frequency spectrum of an ionized sphere. 



The WENO scheme for radiative transfer problems could be incorporated with 
the Euler hydrodynamics. For instance, it is straightforward to generalize the 
calculation of a Stromgren sphere to take account of the hydrodynamic effects 
of the hydrogen gas sphere even when the gaseous sphere undergoes shocks. 
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A Radiative transfer equations 

Considering a universe described by a Robertson- Walker metric, we liave 
ds"^ = g^ydx^dx" = dt^ - a'^{t)gijx'x^ , (A.l) 

wliere we take c = 1. Tlie summation is over to 3 for repeated Greek indices, 
to 3 for Latin indices, a{t) is the cosmic factor and Qij = Sij — hij. We consider 
below the fiat universe, i.e. hij = 0. The distribution function f(t,x\p°') of 
photons is defined by 

dN = -/(x",p")25^(p>„)rfV,rfVp, (A.2) 

where dN is the number of photons in the invariant phase space volume ele- 
ment of 

dVx = a^p^dx^dx^dx^ dVp = dpodpidp2dp3. (A. 3) 

The variables x* is comoving coordinate, p^ = dt/dX, and p* = a{t)dx'^/dX are 
4 momentum, A being the affine parameter. The Dirac-delta function S^{p°'pa) 
in eq.(A.2) is due to the 4 momentum p" of photons to be null, i.e. PaP"" = 0. 

The Boltzmann equation of the distribution function /(x",p") is given by (e.g. 
Bernstein 1988) 

dl^dx^dl_a^,dl^^^ 
dt dt (9x* a (9p* 

where Q is collision term, of which the variables are also t, a;*,z/, n*. Define 
n* = a{dx^/dt), it is a unit vector in the direction of photon propagation, 
n- n = 1. With n*, we have p* = p^n^, or p*^ = hu/c, the frequency u of 
photon. From eq.(A.4), the equation of the distribution function /(t, x^, u, n*) 
is 

where the Hubble constant H = a/ a. From the definition eq.(A.2) of /, the 
specific intensity J oc ly^f, which is also a function of t, x, z/ and n\ Moreover, 
if Q is mainly given by absorption and emission sources, Eq.(A.5) yields the 
equation of specific intensity as follows 
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where k,y is the absorption coefficient, and S the source function. Eq.(A.6) is 
the equation of radiative transfer in the universe. The equation (A. 6) can be 
rewritten as in the format of a conservation flux as 



where a; = In l/u. Therefore, it can be solved by the WENO scheme. 
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